An Euler-Bernoulli-Type Beam Model of the Vocal Folds for Describing Curved and Incomplete Glottal Closure Patterns

Incomplete glottal closure is a laryngeal configuration wherein the glottis is not fully obstructed prior to phonation. It has been linked to inefficient voice production and voice disorders. Various incomplete glottal closure patterns can arise and the mechanisms driving them are not well understood. In this work, we introduce an Euler-Bernoulli composite beam vocal fold (VF) model that produces qualitatively similar incomplete glottal closure patterns as those observed in experimental and high-fidelity numerical studies, thus offering insights in to the potential underlying physical mechanisms. Refined physiological insights are pursued by incorporating the beam model into a VF posturing model that embeds the five intrinsic laryngeal muscles. Analysis of the combined model shows that co-activating the lateral cricoarytenoid (LCA) and interarytenoid (IA) muscles without activating the thyroarytenoid (TA) muscle results in a bowed (convex) VF geometry with closure at the posterior margin only; this is primarily attributed to the reactive moments at the anterior VF margin. This bowed pattern can also arise during VF compression (due to extrinsic laryngeal muscle activation for example), wherein the internal moment induced passively by the TA muscle tissue is the predominant mechanism. On the other hand, activating the TA muscle without incorporating other adductory muscles results in anterior and mid-membranous glottal closure, a concave VF geometry, and a posterior glottal opening driven by internal moments induced by TA muscle activation. In the case of initial full glottal closure, the posterior cricoarytenoid (PCA) muscle activation cancels the adductory effects of the LCA and IA muscles, resulting in a concave VF geometry and posterior glottal opening. Furthermore, certain maneuvers involving co-activation of all adductory muscles result in an hourglass glottal shape due to a reactive moment at the anterior VF margin and moderate internal moment induced by TA muscle activation. These findings have implications regarding potential laryngeal maneuvers in patients with voice disorders involving imbalances or excessive tension in the laryngeal muscles such as muscle tension dysphonia.


Introduction
The configuration of the vocal folds (VFs), a cornerstone of voice production, is determined by the particular combination of activated intrinsic and extrinsic laryngeal muscles. Nominally, the VFs are completely adducted prior to the onset of phonation, and their interaction with the air flow driven by the lungs results in vibrations and consequent acoustic waves, which forms the basis of voiced speech. In some scenarios, complete glottal closure is not attained, which can result in inefficient voice production (Zañartu et al., 2014), and, in some cases, stress concentrations in the VFs that may lead to VF trauma (Dejonckere and Kob, 2009). Hence, incomplete glottal closure is often linked to disorders that are associated with inefficiencies in, or damage to, the vocal mechanism, including Parkinson's disease (Hanson et al., 1984) 1 and muscle tension dysphonia (MTD) (Morrison and Rammage, 1993) 2 .
1 Parkinson's disease is a relatively prevalent disorder, affecting the human central, peripheral, and enteric nervous systems (Braak and Braak, 2000).
2 MTD (Morrison and Rammage, 1993), also known as non-phonotraumatic hyperfunction , is a class of voice disorders associated with misuse of the vocal mechanisms without the presence of organic changes in the vocal organs, leading to low speech quality and vocal fatigue, with a wide range of symptoms and patterns, including excessive/unbalanced activation of intrinsic and extrinsic laryngeal muscles (Roy, 2008; Hocevar-Boltezar Incomplete glottal closure 3 comes in various patterns (Morrison and Rammage, 1993;Nguyen et al., 2009;Södersten et al., 1995) as shown schematically in Figure 1, including bowed shape: a glottal pattern wherein the left and right VF geometries are convex with a gap at the mid-membranous portion; posterior glottal opening: full glottal closure is achieved in the anterior and mid-membranous regions only, leaving the posterior margin open; and hourglass glottal configuration: a pattern with anterior and posterior gaps and potential VF contact in the mid-membranous region. There exist other incomplete glottal closure patterns, sharing similarities with those mentioned above, such as spindle-shaped glottis, and anterior opening (see Rajaei et al. (2014); Södersten et al. (1995)). These latter patterns are not addressed in this work. There exist several experimental and clinical studies in the literature that attempt to elucidate, at least in part, some of the laryngeal mechanisms associated with curved and incomplete glottal closure patterns. Based upon inspection of cadaver larynges, Morrison and Rammage (1993) posited that a posterior glottal opening is associated with excessive activation of the posterior cricoarytenoid (PCA) muscle. Choi et al. (1993b) conducted experimental investigations of excised canine models and found that activating the thyroarytenoid (TA) muscle, while keeping other adductory laryngeal muscles inactive, leads to anterior and mid-membranous glottal closure, whereas the posterior glottis stays open, thus resulting in a closure pattern similar to a posterior glottal opening. On the other hand, they found that when the TA muscle et al., 1998), supraglottal compression (Morrison and Rammage, 1993), and abnormal fundamental frequency (Nguyen et al., 2009;Altman et al., 2005). 3 In this study we refer to the glottal configuration of the VFs at rest immediately prior to phonation initiation, identifying any gaps between the folds as incomplete glottal closure. In clinical settings, incomplete glottal closure typically refers to gaps between the folds when the VFs are at their maximum glottal closure phase during phonation (Södersten et al., 1995;Nguyen et al., 2009). Our definition herein isolates laryngeal factors, which are the focus of this study, by dismissing the dynamics of VF vibrations. is relaxed and other adductory muscles (LCA/IA) are activated, closure is achieved only at the posterior margins of the VFs with mid-membranous opening (Choi et al., 1993b;Chhetri and Neubauer, 2015), thus leading to a bowed configuration as seen in Figure 1. Moreover, complete glottal closure of excised canine larynges was attained via co-activation of all adductory muscles. More recent clinical investigations using refined experimental setups (see, e.g., Chhetri and Neubauer (2015)) further confirm these observations. Chhetri et al. (2012) conducted a parametric study of the effects of intrinsic laryngeal muscle activation, modulated by graded stimulation, on the pre-phonatory posture of a canine model. In addition to confirming the findings of Choi et al. (1993a), Chhetri et al. (2012) found that when keeping the TA muscle activation at a constant level and increasing activation of the cricothyroid muscle (CT), the glottal area increases and the medial bulging caused by the TA muscle activation is reduced. In addition, they observed that when keeping the LCA and IA muscle activation at constant levels and increasing CT muscle activation, the glottis starts to open posteriorly and the glottal area increases. Interestingly, the authors found that with certain muscular executions involving co-activation of the LCA, IA, and TA muscles, the glottis exhibits an hourglass shape (see Chhetri et al. (2012, Fig . 10)).
Besides the aforementioned clinical and experimental works, there exist numerical studies that shed some light onto curved and incomplete glottal closure patterns. Hunter et al. (2004) developed one of the early three-dimensional VF posturing models, where adductory and abductory muscles are incorporated, showing that full activation of the LCA muscle induces nonuniform curvature of the medial surface. Dejonckere and Kob (2009) studied the influence of incomplete glottal closure patterns (with linear and curved VF geometries) on VF vibrations using a multi-mass model, showing that some resting incomplete glottal closure configurations may induce localized VF impact, which they hypothesized to be a potential underlying mechanism inducing VF trauma, especially in females. However, the authors did not study the laryngeal maneuvers that induce these resting glottal shapes. Zhang (2014, 2016) conducted numerical simulations using high-fidelity numerical models, showing that posterior glottal opening occurs with the sole activation of the TA muscle, mid-membranous opening when the LCA and IA muscles are co-activated (without incorporating the TA muscle), and full glottal closure when all adductors are co-activated, in agreement with the aforementioned clinical observations. In a more recent study, Geng et al. (2020) proposed a detailed physiologically accurate finite-element posturing model, based on MRI scan images of a canine larynx. Even though the study does not study the glottal geometry, it provides useful insights into how synergistic activation of laryngeal muscles exhibits complex interaction with laryngeal variables (e.g., VF strain, rotation and translation of arytenoid cartilages, and glottal area). Recently, research interest has also been directed towards investigating how activating laryngeal muscles alters the VF medial surfaces (Pillutla et al., 2022).
Despite these valuable efforts, a clear picture of the physical mechanisms inducing different glottal patterns remains elusive. It is challenging to isolate and control the factors underlying posturing mechanics experimentally. Moreover, high-fidelity numerical models, despite their accuracy in replicating physiological laryngeal postures, do not provide clear intuitive understanding of the mechanics of posturing, and typically suffer from high computational costs. We hypothesize that the non-homogeneous structure of the VFs, which comprise overlapping tissue layers with different mechanical and geometrical properties (Titze and Alipour, 2006), underlies, in part, the different glottal shapes displayed in Figure 1. As such, we propose an Euler-Bernoulli composite beam model of the VFs to elucidate some of the mechanisms underlying glottal patterns prior to phonation. Beam models have been utilized previously to explore VF vibrations and phonation fundamental frequency Zhang et al., 2007). We opt for this relatively simple modeling framework to facilitate exploration of the mechanisms underlying the resulting glottal shapes. To gain refined physiological insights into how intrinsic laryngeal muscles may influence glottal geometry, the proposed model is integrated with the muscle-controlled posturing model of Titze and Hunter (2007).
The organization of this work is as follows: a detailed derivation of the composite beam VF model is introduced in Section 2; analysis is conducted in Section 3; numerical simulations of the integrated beam and posturing model are presented in Section 4; Section 5 presents discussion of the results; and the study is concluded in Section 6.

Model development
Herein, we propose a static Euler-Bernoulli-type composite beam model (see, for example, Bauchau and Craig (2009)) for the VFs, with the different VF layers represented by strata in the beam. For simplicity we assume symmetry with respect to the medial plane; hence, we consider only one (the left) VF. A schematic representation of the composite beam model is shown in Figure 2. The VF beam model consists of three layers: (1) the mucosa with depth d muc and cross-sectional area A muc ; (2) the vocal ligament with depth d lig and cross-sectional area A lig ; and (3) the thyroarytenoid (vocalis) muscle with depth d ta and cross-sectional area A ta . For the sake of compact presentation, we define the index set where muc, lig, and ta refer to the mucosa, ligament, and TA muscle tissue, respectively. We assume that each layer has a uniform rectangular cross-section and the layer thicknesses (in the inferior-superior direction) are equal and denoted by b; thus, layer depth can be computed as Our modelling framework assumes that VF deformation consists of (a) potentially large longitudinal stretching/compression with uniform strain, and (b) modest bending due to the induced moments inside the VF. Let L 0 denote the resting VF length and L denote the VF length after longitudinal deformation due to the associated nominal uniform strainε; that is, We assume the nominal strain is known a priori 4 . Let x ∈ [0, L] denote the position along the deformed VF configuration (after applying strainε) relative to the anterior VF margin, and r denote the depth position along the axis perpendicular to the VF axis relative to the base of the TA muscle (see Figure 2). Consider a plane VF cross-section at position x, and let y muc denote the relative position along the r-axis with respect to the geometrical center of the mucosa (i.e., y muc = 0 corresponds to the geometrical center of the mucosal cross-section). Similarly, let y lig and y ta be analogous coordinates for the ligament and TA muscle, respectively (see Figure 2). Note that the range Let w(x) denote the transverse deflection of the beam (in the r-direction). Moreover, let u i (x, y i ), i ∈ I, denote the longitudinal displacement of the i th VF layer, where longitudinal displacements are with respect to the deformed VF configuration underε. In addition, letū i (x) = u i (x, y i = 0), i ∈ I, denote the longitudinal displacement at the center of the i th layer. Under Euler-Bernoulli beam theory (see, for example, Bauchau and Craig (2009)), the longitudinal displacement functions can be written as where the prime symbol denotes differentiation with respect to x. Continuity of displacement fields necessi- Given longitudinal displacement in the i th layer with respect to the deformed configuration underε, the total strain in that layer is given by 5 Substituting Equation (3) into Equation (5) results in The stress field is estimated from strain and, in the case of the TA muscle layer, TA muscle activation a ta , which is a non-dimensional parameter, ranging between 0 and 1, that corresponds to the activation level in the TA muscle, with 0 indicating a completely flaccid muscle and 1 being maximum contraction. Herein, we utilize local linearization about the nominal strainε. That is, the stress functions in the VF layers, σ i , i ∈ I, are given by the approximate relations where andσ i , i ∈ I denotes the nonlinear stress function associated with the i th layer.
VF tissues exhibit a highly nonlinear hysteretic viscoelatic behaviour (Min et al., 1995;Chan and Titze, 1999). The literature is rich in various studies attempting to develop VF constitutive models that capture, at least in part, the complex mechanical behaviors of the VF tissues (see the review study of Miri (2014) Herein, and for simplicity, we assume the constitutive stress-strain relations associated with the nonlinear stressesσ i to be elastic (functions of strain only) and of exponential type (see Hunter and Titze (2007)) with symmetry about zero strain. In particular, and, in the case of the TA muscle, we include stress induced by muscle activation, resulting in where m i , n i , i ∈ I, are parameters of the constitutive relations, and σ a,max is the maximum active stress Table 1: Numerical values of the geometrical and mechanical properties for each layer in the composite VF model: muc (mucosa), lig (ligament), and ta (thyroarytenoid). Cross-sectional areas are adopted from Titze and Alipour (2006), whereas the parameters m i , n i , and σ a,max are tuned to match experimental stress-strain curves from cadaver and canine models presented in Titze and Alipour (2006, Figure 2.17, p. 88 values of the constitutive relation parameters adopted in this study are listed in Table 1.
The normal forces in the VF layers are computed as Substituting Equations (7) and (6) in yields where denote the nominal normal forces generated by each layer. The total internal normal force is then The moment about the center of the i th layer due to the stress developed in that layer is given by Substituting Equations (7) and (6) into this formula gives where denotes the area moment of inertia of the i th layer.
Let the positions of the geometric centers of the VF layers along the r-axis be denoted r i , i ∈ I; that is, see Figure 2. Take a cross-section at longitudinal position x and consider an arbitrary point on the crosssection located at a vertical position r = r c (see Figure 2). The moment at r c , denoted M c , is given by Consider an element of infinitesimal longitudinal length dx with left edge at position x (see Figure 3), and let V (x) and q(x) denote the shear force and distributed load per unit length, respectively. The force and moment balances on the infinitesimal element yield where second order and higher terms are omitted.
From Equation (14) we deduce that the total normal force N is constant through the VF length. By the assumption that the VF undergoes compression/elongation with associated strainε (see Equation (2)), the force N should be equal to the force that results in that strain, which is the sum of nominal forces F i,0 , i ∈ I, (in Equation (9)). That is,  (8) and (10) into Equation (17) As we are interested in transverse deflection, we aim to obtain a balance equation solely in terms of w.
For convenience, we define and α muc = −l muc , From the continuity condition in Equation (4) and the zero force condition in Equation (18), it can be deduced that the displacement functionsū i , i ∈ I, satisfy the relations Substituting Equation (21) into Equation (13), we obtain where is the composite bending stiffness and is the nominal moment at r = r c due to the nominal normal forces. Note that the bending stiffness is strain-dependent, which can be of importance in posturing scenarios with large VF strains.
Combining Equations (15) and (16) results in which, in terms of the deflection w (obtained by substituting in Equation (22)) is The distributed load q is due to VF contact, which is assumed to be proportional to the transverse overlap beyond the medial plane. That is, where K col is a stiffness coefficient associated with VF contact, θ G is the clockwise angle between the medial plane and the deformed VF configuration under strainε (see Figure 2), and H is the Heaviside function.
In this work, we assume zero transverse deflection at the anterior and posterior ends of the VF. That is, Moreover, we assume zero moment at the posterior VF margin, Furthermore, we assume a reactive moment at the anterior VF margin that is proportional to the rotational displacement with respect to the VF angle at rest, θ 0 . The total angle at the anterior margin between the medial plane and the VF is approximately given by θ G − w ′ (0). Consequently, the moment boundary condition at the anterior margin is given by where K r is a rotational stiffness coefficient. Like the nominal strainε, it is assumed that the angle θ G is known a priori. Finally, we assume that r c corresponds to the geometrical center of the ligament, that is r c = r lig . This assumption, in addition to the boundary condition given in Equation (28), implies that the total normal force N is positioned at the geometrical center of the ligament. This can be deduced from the fact that M c (L) = (r c − r N )N , where r N denotes the r-position of the total normal force N (that is, the force centroid).

Analytical insights from a special case
To gain simple yet useful insights into how internal moments inside the VF beam model affect its curvature, we consider the scenario of zero contact forces (i.e., q(x) = 0) and assume |w ′ (0)| ≪ |θ G |, which reduces the boundary condition given in Equation (29) to Equation (25), with boundary conditions given in Equations (27), (28), and (30), and the definition of M c in Equation (22), can be solved analytically. The curvature of the VF beam model, w ′′ , is given explicitly Note that positive w ′′ implies a convex VF geometry, whereas negative curvature implies a concave geometry.
Recalling the definition of M c,0 given in Equation (24) and implementing the assumption that r By additionally definingM Equation (31) can be rewritten as In the following discussion, we assume that bending stiffness µ c is always positive (µ c according to Equation (23) changes with the elongation or compression of the VF beam model). First, let us analyze abstractly the effects of the moment termsM ta ,M muc , andM r on the VF curvature. We can observe from Equation (35) that the effect of the reactive momentM r on the VF curvature decays linearly with a maximum effect (in magnitude) at x = 0 and zero effect at x = L. In contrast,M ta andM muc have spatially invariant (i.e., constant) effects on w ′′ . The curvature is positively correlated withM r andM ta +M muc .
That is,M r > 0 andM ta +M muc > 0 implies positive curvature (i.e., convex VF geometry), whereasM r < 0 andM ta +M muc < 0 implies negative curvature (i.e., concave VF geometry). Considering the fact that the effect of the anterior reactive momentM r decays linearly along the VF length and the nominal moments induced by the VF layers,M ta andM muc , are spatially invariant, there can arise an interesting scenario for which the curvature changes sign along the VF length. In particular, wheñ w ′′ is positive on [0, x cr ), where and negative for x ∈ (x cr , L], a change from convexity to concavity. The conditions on the internal moments and resulting VF shapes from this analysis are summarized in Table 2. We note that a convex VF geometry is a defining characteristic of the bowed VF pattern. Moreover, the concave VF geometry can be associated with posterior glottal opening. Furthermore, transition along the VF length from convex to concave resembles the hourglass glottal pattern (see Figure 1). This demonstrates that the beam model has the capacity to produce the experimentally-observed glottal configurations shown in Figure 1. Table 2 to physiological posturing scenarios. The termM r , as seen from its definition in Equation (34), is related to VF adduction and abduction, whereinM r > 0 corresponds to VF adduction (θ G < θ 0 ) andM r < 0 corresponds to VF abduction (θ G > θ 0 ). The termM muc +M ta , which is defined according to Equations (32) and (33), is determined by the reactive moments developed in the VF layers, especially the mucosa and TA muscle, during VF tensioning or compression. Note that, based on the area measurements (Table 1) and the assumption of uniform thickness b, the moment arm of the TA muscle, (d ta +d lig )/2, is larger than that of the mucosa, (d muc +d lig )/2. In the case of VF compression, the compressive forces in the TA muscle are typically larger in magnitude than that in the mucosa (i.e., F ta,0 ≪ F muc,0 < 0); hence, the moment induced by the TA muscle,M ta , is positive and predominant makingM ta +M muc > 0.

Now, let us relate the findings listed in
On the other hand, when the VF is tensioned due to activating the TA muscle, the force F ta,0 is positive and predominant and, consequently, the termM ta is negative (see Equation (32)) and predominant. This scenario results inM ta +M muc < 0. These relations between the moment terms and corresponding laryngeal posturing scenarios are summarized in Table 3.
The combined findings presented in Tables 2 and 3 can be summarized by following observations: The bowed shape with convex VF geometry can be due to (a) positive reactive moment at the anterior margin Moreover, the concave VF shape arising in the case of posterior glottal opening can be due to (a) negative reactive moment at the anterior margin (M r < 0) during VF abduction, and/or (b) sufficiently large activation of the TA muscle, whereinM ta +M muc < 0. The hourglass shape may necessitate a coordinated laryngeal maneuver that involves sufficient TA activation (M ta +M muc < 0) and VF adduction (M r > 0) such that M r +M ta +M muc > 0.

Simulations of the combined beam and posturing model
In this section we further investigate VF curvature and incomplete glottal closure by combining our beam model with the VF posturing modeling introduced by Titze and Hunter (2007). In particular, we adopt the implementation of Alzamendi et al. (2022). The posture model relates activation of the five intrinsic muscles to the prephonatory glottal configuration, and in particular, the rotational and linear displacements of the cricothyroid joints and the arytenoid cartilages, where the VFs and intrinsic muscles are modelled as spring-like elements. From the aforementioned displacements, the VF nominal strainε and glottal angle θ G are estimated. Similar to the muscle activation parameter a ta embedded in the VF beam model, the posture model relies on five normalized muscle activation parameters, a ta , a ct , a lca , a ia , and a pca , which correspond to the TA, cricothyroid (CT), lateral cricoarytenoid (LCA), interarytenoid (IA), and PCA muscles, respectively.
In this study, we assume that the muscle activation parameter a ta embedded in the VF beam model is identical in value to the muscle activation parameter a ta in the posturing model (a ta = a ta ).
It is important to mention that the constitutive relations embedded in the VF beam model are different from those in the posture model implementation adopted from Alzamendi et al. (2022). The focus of the current study is to replicate the VF static configurations, wherein we employ experimental stress-strain data   Table 4.
To clearly illustrate the glottal geometries resulting from the simulations, a coordinate system (x 1 , x 2 ) with origin at the anterior margin of the VFs is established. The x 1 -axis is aligned along the medial plane pointing in the posterior direction and the x 2 -axis is perpendicular to the medial plane pointing to the right, relative to the human body frame (see Figure 2). In all figures presented in this section the VF configurations are plotted with respect to this coordinate system; model symmetry is utilized to produce the opposing VF shape.
This section explores several laryngeal maneuvers and how they influence the VF geometry 6 . In par-ticular, and motivated by previous clinical, experimental, and numerical findings (Yin and Zhang, 2014;Morrison and Rammage, 1993;Chhetri and Neubauer, 2015;Yin and Zhang, 2016), we consider laryngeal maneuvers associated with adductory (TA, LCA, and IA) and abductory (PCA), muscles as they have been found to play major roles in inducing curved glottal geometries. The CT muscle has been found to play a major role in regulating phonation fundamental frequency by stretching the VFs, but not in posturing and is thus excluded from this study. Herein, we compare simulation results with findings from previous clinical, experimental, and high-fidelity numerical studies to verify the proposed VF beam model. Moreover, we attempt to elucidate potential mechanisms underlying the curved VF geometries observed clinically by analyzing the beam model details (see Figure 1).
First, we investigate the effects of increasing co-activation of the LCA and IA muscles, which are responsible for adducting the VFs (Alzamendi et al., 2022), while the remaining intrinsic muscles are inactive. Figure 4 presents the glottal shapes and the induced moments corresponding to simulations wherein LCA and IA muscle activation levels are increased simultaneously. Figure 4(left) shows that co-activation of the LCA and IA muscles leads to posterior glottal closure with a remaining mid-membranous gap. This convex VF shape matches previous clinical and numerical findings (Chhetri and Neubauer, 2015;Yin and Zhang, 2016). Figure 4(right) shows that the VF convexity is due to the predominance of the reactive moments at the anterior VF margin (M r is positive and relatively large), which arises due to VF adduction (θ G < θ 0 ), which agrees with the theoretical predictions in Section 3.   Titze and Hunter (2007). Accounting for internal bending moments results in deviation of the VF shape from the linear medial surface prescription (with angle θ G ) of Titze and Hunter (2007). muscle activation levels are increased, while all other intrinsic muscles are inactive. Figure 5(left) shows that isolated activation of the TA muscle leads to anterior and mid-membranous glottal closure with remaining posterior opening, while also shortening the folds. The resulting concave VF shapes are in agreement with previous experimental and numerical investigations (Chhetri and Neubauer, 2015;Yin and Zhang, 2016).
Figure 5(right) shows that the concavity is primarily determined by the internal moments induced by the TA muscle activation (M ta is negative and relatively large in magnitude), which is in alignment with the analysis in Section 3.

Figure 5: (Left) glottal profile and (right) induced moments for increasing TA activation levels and other intrinsic muscles being inactive.
In an effort to explore the mechanics of the hourglass glottal shape, we explore the glottal shape associated with increasing activation of the TA muscle while the LCA and IA are kept at constant non-zero activation levels. Simulating such maneuvers is encouraged by the findings from the theoretical analysis in Section 3 and the experimental observations in Chhetri et al. (2012). Figure 6 displays the glottal shapes and induced moments associated with slight increasing activation of the TA muscle, while the LCA and IA are kept at constant levels (a lca = a ta = 0.6). Figure 6(left) shows that in the case of zero TA activation, the glottal shape is bowed with slight, but not full, posterior adduction. As TA activation is increased, a medial bulge is observed whereas anteriorly the glottal geometry is still convex, resulting in an overall hourglass shape. Figure 6(right) displays how the internal momentsM ta andM muc and the reactive moment at the anterior marginM r satisfy the condition of Equation (36). This aligns with the analysis in Section 3 and suggests that the hourglass glottal shape necessitates involvement of reactive moments at the anterior VF margin (associated with VF adduction) and internal moments induced inside the VF layers (primarily the TA muscle). In addition, this finding is in good agreement with observations in Chhetri et al. (2012), which showed that an hourglass shape is induced by coactivating all the adductory muscles. Figure 6: (Left) glottal shapes, and (right) induced moments for increasing TA activation levels and other intrinsic muscle activation levels being at a lca = a ta = 0.6 and a ct = a pca = 0.
In aggregate, Figures 4-6 indicate that dismissing the TA muscle or the LCA and IA muscles cannot produce full glottal closure, therefore, in the next set of simulations, we investigate the effects of co-activating all of the adductory muscles. Figure 7 shows that (almost) full glottal closure can be attained when all adductors are co-activated (a ia = a lca = 0.45, a ta = 0.7), which is in alignment with previous experimental and numerical investigations (Chhetri and Neubauer, 2015;Yin and Zhang, 2016). Finally, we explore the effects of the PCA muscle, a primary VF abductor, on the glottal geometry, where we consider simulations motivated by the clinical observations highlighted in Morrison and Rammage (1993).

Discussion
The results of Sections 3 and 4 highlight potential mechanisms underlying different patterns of incomplete glottal closure. In particular, results indicate that bowed VF shapes result, in part, from low or null activation of the TA muscle in combination with co-activation of the LCA and IA muscles. The predominant mechanism in this case is the anterior reactive moments that resists bringing the VFs together during adduction. This pattern can also arise in the case of low TA muscle activation and VF compression, as suggested by the analysis in Section 3. In this case, the internal moment induced by the TA muscle tissue (M ta is positive and predominantly large) is the driving factor. This scenario (bowing due to VF compression) can potentially take place when extrinsic laryngeal muscles are excessively activated, especially those associated with VF compressing, such as the thyrohyoid muscle (Hong et al., 1997).
In addition, our analysis suggests that posterior glottal opening with combined VF concavity results from high activation of the TA muscle and low or null activation of the LCA and IA muscles. Our model suggests that the driving mechanism here is the internal moment induced by the TA muscle activationM ta , which is negative and predominantly large in magnitude in this case. A similar glottal pattern also occurs when all adductory muscles are activated in addition to the activation of the PCA muscle. This supports the hypothesis of Morrison and Rammage (1993), regarding the excessive activation of the PCA muscle in patients with MTD. Finally, our analysis suggests that the hourglass glottal shape may emerge from laryngeal maneuvers that involve, for example, moderate co-activation of all adductory muscles, where both anterior reactive moment and internal moment due to TA muscle activation are at play and opposing each other.
The implications above concerning potential connections between incomplete and curved glottal closure patterns and particular muscular executions may help speech therapists to uncover the underlying laryngeal mechanisms associated with some voice disorders. As highlighted in the introduction, incomplete glottal closure can be linked to voice disorders that are characterized by excessive, imbalanced, or deficient activity of the intrinsic and extrinsic muscles such as MTD and Parkinson's disease. Our analysis in the current work suggests two potential mechanisms underlying bowed VFs in some patients with voice disorders (1) the TA muscle is not properly activated (possibly due to muscle activation imbalance), and (2) excessive activation of extrinsic neck muscles, leading to VF compression. Besides, our analysis posits that in patients with abnormal posterior glottal opening and concave VF geometry either (1) insufficiently activate the LCA and IA muscles, whereas the TA muscle is activated sufficiently (in comparison to normal posturing scenarios), or (2) suffer from excessive activation of all adductory and abductory muscles, where the PCA muscle activation mitigates the effects of the LCA and IA muscles, in agreement with the postulation in Morrison and Rammage (1993) concerning patients with MTD. In summary, speech clinicians and therapists may consider the aforementioned candidate underlying mechanisms of curved and incomplete glottal closure patterns when examining patients with voice disorders involving muscular inefficiencies/deficiencies.
A number of simplifying assumptions are embedded in the presented model of this study, including (1) negligible shear deformation, (2) negligible elastic forces from the connective tissues attached to the TA muscle, (3) negligible motions in the superior-inferior direction, (4) neglecting potential bending effects from the vocal ligament by setting r c = r lig , (5) zero moments at the posterior ends of the VFs, (6) small transverse VF deflections, and (7) one-way coupling between the VF beam model and the posturing model, where any contact forces emerging due to the VF curvature do not alter the mechanics of the laryngeal cartilages.
Assumption (1) is a consequence of the adopted Euler-Bernoulli framework. Note that with the uniform thickness assumption, and considering the model dimensions given in Tables 1 and 4, the total VF depth, ∑ i∈I d i , is approximately 10 mm whereas the resting VF length is 15 mm; hence, the depth and length dimensions are quite comparable. For such cases (thick beams), Timoshenko beam theory (Timoshenko, 1921) is typically adopted to account for shear stresses 7 . As the goal of the current work is to construct a simple analytically-tractable model that predicts qualitatively the curved glottal configurations observed 7 It is worth noting that the two theories (Euler-Bernoulli and Timoshenko) do coincide for a uniform homogeneous simply-supported linear beam with specified moments at the end points and zero distributive load, regardless of the beam thickness or mechanical properties (in Section 3, we studied a similar simply-supported case with zero distributive load). These two theories, when compared, tend to produce qualitatively, but not necessarily quantitatively, similar predictions (see, e.g., Beck and da Silva Jr (2011)). clinically, we adhered to the Euler-Bernoulli framework, leaving derivations of more complex models to future work. We posit that assumption (2) is reasonable as the elastic forces from the connective tissues are passive, mostly only restricting the extent to which the VF deflects. Moreover, assumption (3) is suitable as the majority of the VF motion during posturing occurs medially and/or laterally (see, e.g., the findings of Chhetri and Neubauer (2015)). Regarding assumption (4), the ligament is stiffer than other VF layers and it geometrically forms the intermediate VF layer, making it the 'chassis' of the VF layered structure; hence setting r c = r lig , which indicates that the total normal force in the VF is positioned at the center of the ligament (see the end of Section 2), is sensible. Assumptions (5)- (7), in addition to other assumptions such as the rectangular geometries of the VF layers, are introduced to primarily simplify our analysis; hence, further investigation is needed to verify the validity of such assumptions in different posturing scenarios, and refine them when needed.
Despite these simplifying assumptions, our modelling framework is capable of predicting some of the glottal patterns observed in previous clinical and high-fidelity numerical studies, which is encouraging. Still, the speculations and potential explanations provided in this work need further extensive investigation into the biomechanics of VF posturing in both healthy subjects and patients with imbalances or deficiencies in the laryngeal muscles.

Conclusion
In this study, we introduced a simple one-dimensional Euler-Bernouilli-type composite beam model of the vocal folds to understand the mechanisms underlying glottal configuration and incomplete glottal closure.
The model, despite its simplicity, was capable of predicting several clinically observed glottal configurations.
Our analysis highlighted how the different patterns of incomplete glottal closure can arise naturally due to the layered VF structure and the associated induced moments. We coupled the proposed beam model with the posturing model of Titze and Hunter (2007) to gain physiologically relevant insights into the role of laryngeal muscle activation. Our analysis showed that a bowed VF shape can arise due to activation of the LCA and IA muscles without incorporating the TA muscle during adduction or due to VF compression. On the other hand, isolated activation of the TA muscle results in medial bulging and posterior glottal opening.
Posterior opening can also occur due to activating all adductors in addition to activating the PCA muscle.
Moreover, our analysis suggested that an hourglass glottal shape can arise from specific laryngeal maneuvers involving the adductory laryngeal muscles. These results provided potential explanations and conjectures regarding the posturing mechanics of patients with voice disorders such as MTD.
In future efforts we aim to refine our modelling framework, where two-way coupling between the VF beam model presented herein and the posturing model of Titze and Hunter (2007) is incorporated, to account for potential effects that curved VF geometries may exert on the mechanics of laryngeal cartilages. Moreover, we intend to incorporate the beam model with numerical phonation models, to study how curved and partially closed glottal geometries may influence tissue-flow-acoustic interactions, voice quality, and vocal function during phonation.